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We study decoherence in superconducting qubits due to quasiparticle tunneling which is enhanced 
by two known deviations from the equilibrium BCS theory. The first process corresponds to tunnel- 
ing of an already existing quasiparticle across the junction. The quasiparticle density is increased, 
e.g., because of an effective quasiparticle doping of the system. The second process is quasiparticle 
tunneling by breaking of a Cooper pair. This can happen at typical energies of superconducting 
qubits if there is an extended quasiparticle density inside the gap. We calculate the induced energy 
decay and pure dephasing rates in typical qubit designs. Assuming the lowest reported value of 
the non-equilibrium quasiparticle density in Aluminum, we find for the persistent-current flux qubit 
decay times of the order of recent measurements. Using the typical sub-gap density of states in 
Niobium we also reproduce observed decay times in the corresponding Niobium flux qubits. 
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I. INTRODUCTION 

A basic building block of a superconducting quan- 
tum bit is the Josephson junction, which allows coher- 
ent Cooper-pair tunneling. This Josephson current pro- 
vides the necessary non-linear element in an electric cir- 
cuit and enables reduction of the externally controllable 
quantum dynamics to involve only two eigenstatesi. A 
drawback for superconducting qubits has been that these 
systems can be very sensitive to various effects in their 
nearby environments. The influence of two major deco- 
herence sources, charge fluctuations and two-level fluc- 
tuators, have been substantially reduced in recent years. 
Qubits that are less sensitive to changes in background 
charge have been build^-tS^ while fluctuators have been re- 
moved by decreasing the junction size. Recently, trans- 
mons and persistent-current flux (p-flux) qubits with ex- 
tremely long decay times have been demonstratedi*5- . It 
is therefore of great interest to investigate mechanisms of 
decoherence that were previously unobservable. Temper- 
ature dependent decay time measurements in phased and 
transmon qubits^ suggest that non-equilibrium quasipar- 
ticles might have now become the main factor limiting 
the decay times of superconducting qubits. 

In this work we analyze qubit decay and dephasing due 
to quasiparticle tunneling. As the source of tunneling we 
consider the two following deviations from the equilib- 
rium BCS state of superconductors. The first one is a 
non-equilibrium distribution of quasiparticles, which cor- 
responds to a finite quasiparticle density above the gap, 
present at all temperatures^"—. Non-equilibrium quasi- 
particles could originate in a quasiparticle diffusion from 
higher temperature regions in an experiment, or through 
stray radiation, and are relatively long lived due to slow 
quasiparticle recombination rate — "— . The second pro- 
cess we consider is quasiparticle tunneling by breaking of 
a Cooper pair, which can occur at typical energies of su- 
perconducting qubits if there is an extended quasiparticle 
density inside the gap. It has long been observed that the 
number of states inside the gap is vastly larger than the 



predictions of BCS theory^l. This effect is especially pro- 
nounced in Niobium, a material that is widely used in the 
production of qubitai^. Nonetheless, Aluminum based 
superconducting devices are observed to posses sub-gap 
states tooiSiiSj but with considerably smaller density as 
Niobium. 

Our results are that the two quasiparticle processes 
result in a similar type of qubit decay rate. This pro- 
vides a clear comparison between the magnitudes of de- 
coherence due to the two sources. For conventional Alu- 
minum qubits non-equilibrium quasiparticles should have 
a more significant impact on decoherence, for Niobium 
qubits the sub-gap density of states. We also find that 
the sensitivity to quasiparticle tunneling is vastly differ- 
ent for differing qubit designs. By using the lowest re- 
ported density for non-equilibrium quasiparticles in Alu- 
minumi^ we obtain p-flux qubit decay times similar to 
recent experimentSi^, whereas for other qubit types the 
results are from two to three orders of magnitude greater 
than the observed ones. For Niobium p-flux qubits we 
also reproduce usual experimental decay times by as- 
suming a Dynes-type of quasiparticle sub-gap states, ob- 
served in various other experiments^^. These central re- 



Qubit type 



Transmon— 
Phase*^ 

Capacitively-shu. flux^- 
Persistent-current flu»2^ 



Non. QP (Al) 



1.2 ms 
0.6 ms 
80 ^s 
12 ^s 



Sup-gap DOS (Nb) 



20 /is 
14 ^s 
3 /is 
0.3 /is 



TABLE I. Comparison of decay times obtained by assum- 
ing a non-equihbrium quasiparticle density Wqp — 0.033//im"^ 
(second column. Aluminum) or a Dynes parameter 7Nb = 
Fd/A — 10~^ (third column. Niobium). We used qubit pa- 
rameters specified in each of the citations at each qubit. We 
also assume Aai = ANb/7 = 180 /leV and the validity of 
the Ambegaokar-Baratoff relation. Typical sub-gap density 
of states in Aluminum (7A1 ~ 10"'* — 10~^) does not limit the 
coherence here. 
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suits are summarized in Table HI 

The paper is organized as follows. In section |lT] we 
introduce the Hamiltonian and methods describing the 
qubit-quasiparticle dynamics. The section lllll is devoted 
to discussion of decoherence processes and the corre- 
sponding rate equations. We also discuss decay rates for 
single-junction qubits obtained by using typical experi- 
mental parameters of the quasiparticle environments. In 
section ITVl we generalize the treatment to multi-junction 
qubits, where a special concentration is given to flux 
qubits in their different designs. In section |V] we sum- 
marize the results. 



II. FOUNDATIONS 

Our starting point is the division of the total Hamilto- 
nian into three parts, 



total 



qubit 



qp 



(1) 



Here the qubit Hamiltonian TJqubit describes the collec- 
tive degrees of freedom that constitute our two-level sys- 
tem. The quasiparticle Hamiltonian Hqp describes the 
electronic degrees of freedom in the superconductors, 
treated independently of the collective part. Interaction 
between these two parts emerges due to single-electron 
tunneling modeled by the tunneling Hamiltonian Ht. It 
describes, e.g., quasiparticle tunneling across the Joseph- 
son junction which can cause transitions between the 
eigenstates of the qubit. In the following we discuss in 
more detail the central properties of each of the three 
parts. 



B. Quasiparticle degrees of freedom 



The quasiparticle degrees of freedom in superconduct- 
ing leads are modeled by the Hamiltonian Hqp. For all 
further calculations the quasiparticle properties are as- 
sumed to be independent of the qubit state. We consider 
here two models as two separated causes of quasiparticle 
tunneling at energies well below the superconducting gap 
A. We label these as model / and model //. 

Signatures of the electronic structure in the two models 
are found from correlation functions of type, 

GL(k,i) = (Cka(i)4 \ , G<Jk,0 = (cLckaW) , 



(3) 

Here is an electron annihilation (creation) operator of 
the state k with spin a {a (3). The Fourier transforms 
of the correlation functions are related to spectral and 
distribution functions as^l 

GL(k,^) = A(k,^)/±H, 

(4) 

F.g(k,c.) = -F^^{k,L,) - B(k.w)/±H. 

Here ^ f and /+ = 1 — /. In equilibrium / — 
/oq(cj) = -I- exp{ijj/kBT)]. We define the normalized 
density of states n and the pair density p as 



(5) 



A. Single-junction qubit Hamiltonian 

As the general single-junction qubit Hamiltonian we 
use the following sum of capacitive, inductive, and 
Josephson coupling energy 



qubit 



EcN^ + V'ext)' - Ej cos^. (2) 



Here D is the density of states nearby the Fermi surface, 
including spin. These functions, added with the distri- 
bution function /, have their own characteristic forms in 
the two models and are all we need to know about the 
quasiparticle environments in the final forms of the qubit 
decoherence rates (see section IIIII) . 



1. Model I : N on- equilibrium quasiparticles 



Here N is the number operator of electron charges on 
the junction capacitor C. It is a conjugated variable 
to the phase difference across the junction (p, satisfying 
[if, N/2] = i. This choice of qubit Hamiltonian corre- 
sponds to the phase and dc-flux qubits, i.e., when the 
junction is placed in a superconducting loop with induc- 
tive energy El and thread by flux ipcxt- For the trans- 
mon one has El — and periodic commutation relations 
[iV/2, 6=^*"^] = ±e^''^, related to the 2e-quantization of 
the charge that tunnels across the Josephson junction. 
All qubits we consider work in the limit where the Joseph- 
son coupling energy Ej is larger than the single-electron 
charging energy Ec = e^/2C. Generalization to multi- 
junction qubits is straightforward, see section ITVl 



In model / we assume that in each lead quasiparticles 
can be described by the BCS Hamiltonian but exist with 
a general (non-equilibrium) probability /. In the BCS 
model the quasiparticle and pair densities have the form 

n(^) = ^=^6(0.2 - A2)sgn(^), 

A 

piLo) = —==e{uj^ - A2)sgn(c.). 
VtJ^ — A^ 

An important quantity is the total density of quasiparti- 
cles 



(7) 



riqp — 2D / n{uj) J {uj)duj . 
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In calculations n^p is assumed to be a given constant. 
This model lias also been investigated in Refs.[Tlland[T2l. 



2. Model II: Sub-gap density of states 

As the other source of quasiparticles we consider the 
presence of quasiparticle states below the gap. We as- 
sume that the densities can be expressed in the same form 
as for the BCS theory, but with the following broadening, 




V(c^ + «rD)2-(A + zA2)2 

A + zAa 



sgn(w), 



sgn(w). 



(8) 



Here Fd <C A is the so-called Dynes parameter— and 
A2 ^ A is a possible imaginary part of the supercon- 
ducting gap^''. We consider Fd and A2 to be given pa- 
rameters and independent of energy. The effect of these 
modifications of the BCS theory are very similar nearby 
the gap, but differ at low energies a; w (see section 
IIII A2[) . The broadening through an imaginary part of 
the superconducting gap^, A2 (Fd = 0), is supported by 
the microscopic Eliashberg theory. However, a constant 
imaginary part of the energy, Fd (A2 — 0), usually fits 
better to experiments done with, e.g., Niobiu m^°i^^i^^ . 
For all further purposes we will only consider finite Fd 
and A2 = 0. 



C. Quasiparticle tunneling 

The interaction between the qubit and the quasiparti- 
cle environments occurs due to the possibility of quasi- 
particles to tunnel across the junction(s) and simultane- 
ously change the charge configuration in the qubit space. 
This is described by the modified tunneling Hamiltonian 

Ht = Ht+ Ei cos if, (9) 
where the hare tunneling Hamiltonian has the form 



i^T =t^(4aClaT + h.C.). 



(10) 



klo 



Here the states k and 1 belong to the opposite sides 
of the junction. We assume a constant tunneling ma- 
trix element t, related to the normal state tunnel resis- 
tance as i?T = h/t^D^ne^. The charge-transfer opera- 
tor f = J2n \^ + = e**^/^ accounts for the cor- 
responding changes in the charge number of the leads. 
Cooper-pair tunneling, which is a second-order process in 
quasiparticle tunneling, is already included in the qubit 
Hamiltonian in an approximative way as the operator 
—Ej cosif. Therefore we need to subtract it from H^^ to 
avoid double counting. 



III. DECOHERENCE RATES 



The second-order expansion of the qubit's density- 
matrix equation of motion in the tunneling Hamiltonian 
leads to contributions describing qubit decoherence and 
parameter renormalization. Here we discuss in detail the 
decoherence terms, i.e., terms leading to energy decay 
and pure dephasing. We then estimate their magnitude 
for usual single-junction qubits. An exact derivation of 
the rates, a formulation for the energy-level renormaliza- 
tion effects, and a discussion for the effect of higher-order 
terms is given in the Appendix. 

The system we consider is a two-level system interact- 
ing with a fermionic bath. The coupling to the bath can 
be divided into two parts: (i) the part causing transi- 
tions, oc ffx, and (ii) the part causing pure dephasing, 
cx . The decay rate, Fi, is a result of coupling to a^- 
This leads also to the ordinary dephasing rate F2 = Fi/2. 
Fluctuations through cr^ lead to the pure dephasing rate 
F2» , which can be interpret to be a result of low-frequency 
fluctuations in the qubit energy splitting due to coupling 
to the batbi. 



A. Tunneling processes 



In leading order we obtain two distinct types of tun- 
neling processes causing decoherence. The first one is 
tunneling of an existing quasiparticle. In the case of the 
BCS state this is described by an operator of the type 



Oi — MkWiT — VkViT^ . 



(11) 



This process contributes mainly in model /. The second 
one is the breaking of a Cooper pair, described by the 
BCS operator of the type 



O2 — UkVlf + VkUlf^. 



(12) 



This process contributes only in model //. Impor- 
tant here is that the operators Oi are superpositions of 
two electron-tunneling directions, T and T^, as long as 
M, u 7^ 0. As the qubit states are superpositions of differ- 
ent charge states, the direction of electron tunneling be- 
comes indistinguishable and the two tunneling processes 
interfere. This gives rise to a phase dependence in quasi- 
particle tunneling, similar to the cos t^-term in the classi- 
cal Josephson effect^^. This interference effect has also an 
important role in decoherence of superconducting qubits, 
as discussed below. 
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1. Qubit decay 

For the qubits considered in this work the qubit decay 
rate can be written in the form 



oo poo 



oo •/ — oo 



(tlTli) / du duj'f-{uj)f+{uj') 



X [n{uj)n(Lu') + p(lu)p{ll>') cosipo] S{uj — uj' ~ SE). 



(13) 



Here SE is the energy-level splitting of the qubit and the 
angle (po is defined as 



e^" = — - 



(t|ru)a|f| t) 



(t \f\ i) 



(14) 



The matrix elements depend on the choice of supercon- 
ducting qubit, but as we will discuss, they are in fact 
very similar for broad classes of qubits. 

Expressions such as (IT3)) are often derived in the BCS 
excitation picture, where integration over only positive 
(excitation) energies emerges. In this semiconductor- 
type presentation negative energies appear as well. But 
the result of the two presentations is the same. The pos- 
itive integration region of both frequencies a; and w' cor- 
responds to tunneling of an existing quasiparticle to one 
direction, as the negative integration region to the other 
region. The contribution with positive uj but negative 
w' describes breaking of a Cooper pair during the single- 
electron tunneling. It can exist only if there is an ex- 
tended states below the gap (because SE < 2 A). The 
contribution with positive u' but negative lu contributes 
only if SE < and would correspond to recombination 
of two quasiparticles with excitation of the qubit. 



2. Qubit decay m the low-energy approximation 

To simplify equation we assume SE <C A and 
that the width of the quasiparticle distribution (^ ksT) 
is much smaller than the qubit splitting, SE. This means 
that for model / the distribution / is nonzero above the 
gap (w > A) practically only in a very narrow region 
[we have always /(— w) = 1 — /(i^) = /+(a;)]. On the 
other hand, for the model // with Dynes broadening (and 
A2 = 0) one can approximate 



P = 



A ' 
A a"' 



(15) 



For a finite imaginary part A2 and Fp = one would 
have n = (tj/A)(A2/A) and p = {uj / A)^ {A2 / A) . 
However, in the following we assume a Dynes-type 
broadening22i2^i^ when considering model II. 



In the discussed low-energy approximation we obtain 
the following common form for the decay rates. 



(1 + ei cos ipo) , i = I,II . 



(t \T\ i) 



(16) 



Important here is that the rates are proportional to ^i, 
defined for the two cases as 



6 = 



i(A + SE) 



D 



(17) 



^n^SEi^-^ 



The form Fi cx 1 -I- cos ipo for the decay rate (|T6l) is 
a result of the interference effect in quasiparticle tunnel- 
ing2^ (see the discussion in section UlI A[) . Its magnitude 
and sign are given by 



1 



1 + ^£;/A 



SE 
~A 



(18) 



For model / we have e/ « 1 whereas for model // e// w 0. 
Physical interpretation for this is that the subgap states 
in model // are close to metallic states and do not show 
significant pair correlations (p « 0), needed for the effect. 



3. Pure dephasing 

In the leading-order expansion of the qubit's time evo- 
lution the pure dephasing appears as an extra decay of co- 
herent oscillations through transition terms proportional 
to (t |r(t)| t) or a |r(t)| I). The pure dephasing rate 
is then similar to the rate by setting SE — 0, and 
changing the matrix element, 



Fo. - 



dLo'r{^)f^{^') 

X [n{uj)n{uj') -\- p{uj)p{u;') cosf^o*] (5(a; — cj') , (19) 

= (t |f 1 1) - a |T| ;) ' . 



Here we have assumed that 



( a 1^1 i) ] 


2 


f (t |f 1 t) ^ 


\ 


a 1^1 ;) 


) 




(t \f\ t) 





(20) 



which is the case for most of the qubits considered in this 
paper. The only exception is the flux qubit away from its 
symmetry point, discussed more detailed in section ITVl 

The value of the integration in Eq. ([T9| depends on 
the distribution /, which depends on the model used. 
In model // we assume the equilibrium distribution 
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/eq = 1/[1 + exp(a;/fcBr)]. Then the pure dephasing 
corresponds to fluctuations originating from tunneling of 
thermalized (sub-gap) quasiparticles. For this case the 
integration can be easily evaluated giving (r-£,/A)'^kBT. 
This means that for small temperatures ksT <^ SE the 
pure dephasing rate is much smaller than the decay rate. 

In model / the situation is more difficult due to the log- 
arithmic divergence of the energy integral at the energy 
gap. Similar singularities appear also in other properties 
of superconductors^^,, but stay finite due to finite lifetime 
effects or gap anisotropy. It is now crucial, which type of 
broadening is assumed. If one introduces a small imag- 
inary part A2 and accordingly redefines Tiqp in Eq. (O, 
then the integral in Eq. (fT9|) remains small. As a result 
of this the pure dephasing rate stays small compared to 
the energy decay rate. Other methods for circumventing 
the divergence can produce larger results for the inte- 
gral. However, as we will discuss in the next section, 
pure dephasing is additionally suppressed by small tun- 
neling matrix elements. 



B. Discussion for single-junction qubits 

We will now discuss the qubit part of the decoher- 
ence rates in the case of single-junction qubits. Specifi- 
cally we consider the phase qubit, the transmon, the dc- 
flux qubit, the strongly anharmonic phase qubit^ and 
the fluxonium. Strictly speaking the fluxonium is not 
a single junction qubit, but we consider the integrated 
Josephson junction array simply as an inductance. For 
a detailed discussion of the Fluxonium including quasi- 
particle tunneling in the junction array see Ref. Il2l The 
discussion can be structured along two major properties. 
One is the squared magnitude of the qubit matrix el- 
ement in Eq. (|16p or of the difference between matrix 
elements in Eq. (|19p , divided by the tunneling resistance 
i?T- The second is the phase difference (po in the inter- 
ference (cost^o) term. 

When it comes to estimating the matrix elements 
qubits fall broadly in two classes. One class of qubits 
have eigenstates located in a single local minimum of the 
qubit 's potential energy, 



U : 



Er 



(if) - </7cxt)^ - Ei COS ip. 



(21) 



The lowest eigenstates are then similar to that of a har- 
monic oscillator and they are symmetric or antisymmet- 
ric around the local minimum. Here the decay and de- 
phasing elements are in a good approximation given by 



Ec 
SE' 
E^ 
SE 



(22) 



In this class we have clearly the transmon and the phase 
qubit. A qubit with a quartic potential as proposed in 



Ref. |30 (strongly anharmonic phase qubit) has the same 
matrix element and despite its rather unusual potential, 
the fluxonium belongs to this group as well for ipcxt — 0. 
The qubit that bucks the trend is the dc-flux qubit, which 
in fact is very similar to the fluxonium for ipcxt = With 
its double well potential it is very different from all the 
other qubits. We will discuss the form of the matrix ele- 
ments in detail in the next section, but at the symmetry 
point they can be shown to have the approximative form 



71 ,r2 ■ 2 f "if mill \ 

M, = sm j , 



M4 



(23) 



where ipmin is the solution to the transcended equation 

ELPmin - Ej sin <^min = , (24) 

and 



I El - Ejcos{ipr, 



2Ec 



(25) 



The decay matrix element is the largest of the discussed 
qubits. Important is also that this qubit has relatively 
low tunneling resistance i?T. Depending on the interfence 
effect it can be rather sensitive to quasiparticle tunneling. 

The qubits have various interference angles, as defined 
in Eq. (ITi|) . For the qubits that are similar to a har- 
monic oscillator this angle is given by the position of the 
local minimum. The potentials of the transmon and flux- 
onium qubits are symmetric around the phase ipo — 0, 
whereas the potential of the phase qubit is also practi- 
cally symmetric, for certain ipo ^ 0. In the dc-flux qubit 
(with (poxt = tt) the interference angle is given by (fo = tt, 
and the situation is the same for the stronlgy anharmonic 
phase qubit. This provides protection against qubit de- 
cay due to non-equilibrium quasiparticles, as then the 
interference between two electron tunneling directions is 
destructive. However, this protection is only partial be- 
cause for typical qubit splittings SE « A/ 10 one has 
ei w 0.9 (< 1). 

The common decay rate (jl6p shows that the two quasi- 
particle sources produce a similar type of decay rate. It 
is now easy to estimate which source should be domi- 
nant in typical experimental conditions. A crucial quan- 
tity is the dimensionless parameter £,i/A, defined in 
Eq. pT|) . The division by the energy gap also accounts 
for changes in Eq. ^6]) due to change in the tunneling 
resistance, if the same qubit (with same Ej) is build 
from, for example. Niobium instead of Aluminum, as 
one has Rt oc A. For an Aluminum superconductor 
the parameter is usually measured to have a value in 
the range ^//A ^ 10^^ — 10^^ (non-equilibrium quasi- 
particles) or C///A ^ 10^^^ — 10~^ (sub-gap density of 
states). This means that for Aluminum non-equilibrium 
quasiparticles have usually an impact that is several or- 
ders of magnitudes larger than that of sup-gap states. 
On the other hand, for Niobium using the typically ob- 
served Dynes parameter^Si?! we obtain values ^///A ~ 
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FIG. 1. (a) The original persistent-current flux qubit as pro- 
posed in Ref. [3ll . The p-flux qubit has three junctions, two 
large junctions with Josephson energy Ej and phase differ- 
ences ifi and ^3. The smaller junctions has the Josephson 
energy Ejs and the phase difference (fi2- Through the loop 
a flux ify is applied, (b) For technical reasons many flux 
qubits actually have four junctions. The fourth junction has 
the Josephson energy i5,j4 and the phase difference (^34. (c) 
To allow for more tunability often the smaller junction is re- 
place by a SQUID loop consisting of two junctions with phase 
differences ip2a, and ip2b 



large junctions. In another version of the flux qubit, the 
small junction has been replaced with a SQUID, meaning 
by a loop with two junctions. Fig. [TJc). This allows for 
more tunability. 

We will start this section with a discussion of the gen- 
eralization of the decoherence rates to many junctions. 
Then we will continue with calculating the rates for the 
p-flux qubit with three junctions and then later discuss 
the modification by the fourth junction and the tunable 
third junction. We end the section with the capacitively 
shunted flux qubit, which by its properties is shown to 
be similar to an anharmonic oscillator. 



A. Generalization of the decoherence rates 

As before we define our total Hamiltonian of the qubit 
and the quasiparticle environment as 



10~^ — 10~^, corresponding to the upper bound values 
for non-equilibrium quasiparticles in Aluminum. 

In table U we compare decay times for various qubits 
calculated by assuming the lowest reported value for the 
density of non-equilibrium quasiparticles in Aluminum, 
given in Ref. [lol For single-junction qubits (transmon 
and phase) we find decay times of the order of ms, which 
is two to three orders of magnitude higher than the best 
experimental observations. For the same qubits but build 
from Niobium the decay rate, Eq. (I16p . simplifies to 



Fi 



1 



CRt{A/Tu)^- 



(26) 



being consistent with the result of the classical limit. In 
this case we find decay times of the order of 10 fis. 

The single-junction qubit that is the most sensitive to 
both quasiparticle sources is the dc-flux qubit, for which 
we obtain decay times of order 100 fjs (Al) and 100 ns 
(Nb) (not listed in table HI). The value for Aluminum 
includes the protection due to destructive interference 
which is approximately a factor of 10. All together, 
in these conditions the Aluminum based single-junction 
qubits are protected against non-equilibrium quasiparti- 
cles, most of them due to small i?c/i?j-ratio. 



IV. MULTI-JUNCTION FLUX QUBITS 

In this section we generalize our analysis to multi- 
junction qubits. We consider in detail the most common 
of them, the fiux qubit, which comes in many different 
shapes and forms. The original proposal. Fig. [TJa), is a 
qubit with three junctions^, but in general it is today 
used in many variations often with several extra junc- 
tions. For technical reasons the flux qubits have always 
been build with a fourth junction. Fig. [Ijb), which in 
most modern flux qubits is of the same size as the two 



qubit 



qp 



(27) 



where iJqubit is the Hamiltonian of the flux qubit, cou- 
pled via the tunneling Hamiltonian to a quasipar- 
ticle environment Hqp. With a total of n junctions we 



have Ht = Ht - 
Hamiltonian is 



up J 



where the bare tunnel 



H 



(28) 



klo 



Here c^j^'' and c^j'' are electron annihilation (creation) op- 
erator of the states k and 1 on the two sides of the junction 
j with Josephson coupling Ejj , and Tj is the correspond- 
ing qubit-space charge-transfer operator, given in differ- 
ent cases below. The quasiparticle environments in each 
of the leads are assumed to be identical and described by 
relations PI5)) . 

In the leading order the decay rates can be written into 
a similar form as for the single-junction case, Eq. (jl3p . 
and in the low-energy approximation fsection llll A2[) one 
gets 



"pt( 
^ 1 



total 



2^2 



(1 + Ci COSCpoj) , 



(29) 



Ml = 



(t|T,U) 



Here Rtj and Lp^j stand for the tunneling resistance and 
the interference angle of the junction j. Similar general- 
ization applies also for the pure dephasing rate, Eq. pj|) . 
In the following we discuss the form of the matrix el- 
ements Mfj and the phases ip^j in the considered flux 
qubit realizations. 
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B. The persistent-current flux qubit 

The first persistent-current flux qubit we consider has 
three Josephson junctions, see Fig.[TJa). Its Hamiltonian 
can be presented in the form 



qubit 



Er 



I 



-N'i + Nt 



1 + 2as 

- Ej [2 cos((p+)cos(cp_) + acos(2<y9+ - (pcxt)] ■ 

(30) 

Here Ec is the total charging energy of the system, 
Ej is the Josephson energy of the larger junctions and 
a = Ejs/Ej is a parameter that determines the energy 
splitting of the qubit. In the following it is assumed that 
for the corresponding capacitive term it holds = a. 
The ratio of Josephson and charging energy for the p- 
flux qubit is generally given by Ej/Ec ~ 50. We have 
three junctions with phase differences ipi, ip2 and 1^3, as 
illustrated in Fig. [ija). The phases in the Hamiltonian 
are then defined by 



1 



(Vl ± Vs) , tPoxt = Vy 



(31) 



For each dynamical phase we have a conjugate charge 
variable, [<^±,-/V±/2] = i, and for each junction we have 
a tunneling operator given by 



T3 = g''('P+-¥'-)/2 ^ 

1 . Eigenstates 

The potential of the p-flux qubit is given by 

U{ip-,ip+)/ Ej = — 2 cos((^_|_) cos((^_) 
- acos(2v3+ - (^oxt)- 



(32) 



(33) 



In the direction of the system behaves similar to a 
harmonic oscillator that remains in the ground state. 
Therefore we will focus on the dynamics in the ipj^ di- 
rection. Choosing (^cxt = (the symmetry point) and 
(y9_ = we expand U for small (y9+, 



U/Ej « (1 - 2a)ifl + 



2a 

T 



1 

12 



(34) 



We have a > 1/2. This means that the eigenstates of 
the qubit are formed by symmetric and antisymmetric 
superpositions of the ground state of two wells. These 
ground states are centered around the minima ±(^+ min, 



^/^^/2c 



1 



V8a- 1 

In the left well the eigenstate is of the form 



(35) 



1/4 



= - exp [-c.'((^+ + ¥'+,min)V2] , (36) 



where we defined 



(4a2 + 1) 



El 

2E, 



c 



(37) 



The state in the right well, has the same form, just 
centered at (^s+.jnin- Under the same approximation we 
used to derive the potential given by Eq. ([M)) we can 
now find the coupling strength t between the two states, 
which is given by 



- = (i?|7Jqubit|i) 

[-2 - 3a + 2a2 - 2^1 ^^{^a" + 1)]^,; 



(38) 



1 + 2a 

X exp(-wV+,min)- 

Here we assumed ipc^t = tt and consider the contribution 
of the order of Ej . There is an additional contribution of 
the order of ^/EjEc and Ec, but they remain small. Of- 
ten the coupling strength has been calculated using the 
WKB approximation. Our calculations underestimates 
the coupling strength but still gives a good order of mag- 
nitude approximation. 

We move away from symmetry like ipcxt = tt — Sipcxt 
and in the lowest order of Sipcxt we get 



^qubit ~ -ff qubit ("^cxt = '^)+Ej ( -~2aip 



Aaip^ 



(39) 

If we now compare the energy of the minima with the 
energy of the minima for (5(/?ext = we find the energy 
difference 



^ ^ 12aV6V2^^, 

(8a -1)3/2 '^^ext. 



(40) 



Now we can write the Hamiltonian in the standard way 
of a two-state approximation for the states \L) and \R), 



I 5e t 



qubit 



2 y t -5e J 
The basic form of the states is therefore given by 



t)=cos(-l|L)+sin(-)|i?), 



(41) 



(42) 



U> = -sin(^-j |L)+cos(^-J \R), 

with tan 9 — t/Se. 

Let us now consider an operator of the form 



T 



We get in a good approximation 



{L\f\L) 



-jV+,mi„/2 



{R\f\R) 



and 



{L\f\R) - {R\f\L) = e- 



(43) 

■+.™. (45) 



Using these equations we will be able to calculate the 
effect of quasiparticle tunneling on the fiux qubits. 
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2. Qubit decay 



3. Dephasing 



To estimate the qubit decay rate we start from the 
transition element of the operator T = e'^'^+^^. This can 
be expressed in the form 



The dephasing is characterized by the diagonal matrix 



a 1^1 1> = 



I sm a sm 



+ cos 9e 



-.(46) 



From this general form we can deduce the results for the 
tunneling operators Ti and T3, 



(;|Ti|t) = a|T3|t) =*sin(?sin 
+ cos6le"'^''^+.-», 



(47) 



where we have used the fact that the system stays in the 
ground state in the direction (^_. Important here is that 
for the symmetry point one has = 7r/2 and the matrix 
elements (|47p are purely imaginary. This is in contrast 
to the (single-junction) dc-flux qubit where the element 
is real. It follows that the interference angles in Eq. (|29|) 
are (poi — <^03 = 0, i.e., one has constructive interference 
instead of destructive. 

The junction that is somewhat different is the second 
junction, for which the matrix element is given by 



{i\T2\^) = 



— 2sm6'sin(^-) 



cos ae 



gi'i'oxt/2 

(48) 

The element (I48p has an additional phase factor, g*"^""*'^, 
and becomes real at the symmetry point (/Scxt — This 
corresponds to the case of the dc-flux qubit and destruc- 
tive interference. 

The interference effect has an important role in de- 
termining the sensibility to non-equilibrium quasiparticle 
tunneling. As discussed in section llll Bl the decay matrix 
element of this qubit type is the largest, but the single 
junction dc-flux qubit is protected by destructive interfer- 
ence. As a result of constructive interference for two junc- 
tions of the p-flux qubit, this protection is lost. It follows 
that even with the smallest observed non-equilibrium 
quasiparticle density^, corresponding to ^//A ^ 10~^, 
we obtain a decay time of the order 10 /iS. This result 
coincides with recent experiments^. 

Results for Niobium p-flux qubits is similar to dc-flux 
qubits: using typical subgap densitiesSS. one obtains de- 
cay times of the order 100 ns, being also similar to recent 
experiments^. These results are compared with the ones 
of other qubit types in table ID 



elements of T., , which can be shown to have the form 



(t iTil t> = (t 1^3 1 t> - -zcos0sin(^^±p) 



sm( 



°+,mi„ + COS 



a |Ti| ;) = a IT3I ;) = +*cos0sin [^^) 



sin6'e"" '^+-""° + cos 



(49) 



(50) 



The second junction has again an extra phase depen- 
dence, 

±z COS 9 sin mm) =t sin 9e ^ '^+,niin _^ [(p^^^:^^) 

(51) 

The presence of a term that does not change sign when 
changing the state leads effectively to a different angle 
(/?o* for the two states | f) and | J,), and deviations from 
result (19). However, near the symmetry point such cor- 
rections stay small and one can neglect this effect. In this 
region the main contribution to pure dephasing comes 
from the junction 2. Using Eq. (ICTI) one obtains then for 
the relevant dephasing matrix element 



|(t|r2|t)-a|T2U)|' = e-2'^''^+- 



(52) 



This is usually much smaller than the decay element. 
Combined with the discussion of section llll Bl this means 
that at the symmetry point the dephasing is limited only 
by the qubit decay. 



C. The four-junction flux Qubit 

For technical reasons the p-flux qubit is in fact always 
build with four junctions instead of the necessary three. 
The fourth junction has a charging energy Eca and a 
Josephson energy iJj4. We can assume i?j4 ^ Eci- The 
additional junction adds another dimension to the prob- 
lem, and we can approximate states along this dimension 
as states of a harmonic oscillator. Within this approxi- 
mation the total Hamiltonian can be written as 

i?qubit = \5Et^ + g (miTj; -I- m2T^) {a^ + a) 

2, (5d) 

-I- wa^a, 

where Ti are the Pauli matrices acting on the eigenstates 
of the p-flux qubit ([32]), with the energy splitting 5E = 
^/5e^ + 1'^. These states are coupled with the coupling 
strength g = Eja{2EcA/ Ej4)'^/^ to the fourth junction, 
modeled as a harmonic oscillator with frequency w = 



9 



y/WEciEjl. The couphng is determined by the foUowing 
matrix elements, 

mi = (t I sin {2ip+ - (^cxt) | 1), 
2m2 = (t I sin (2(/3+ - (^ext) I t> - (I I sin (2(p+ - ip^^t) \ \.) 

(54) 

If the fourth junction is large, £'j4 — > oo, the couphng 
goes to zero. In this case the fourth junction is decoupled 
and always stays in the ground state. However, in many 
flux qubits the fourth junction has the Josephson energy 
i?j4 ~ In this case the couphng is strong, g k uj. 

The tunnehng operator for the fourth junction is sim- 
ply given by T4 = e^Vil"^ ^ Close to the symmetry point 
we have mi k, and additionally we assume w 3> In 
this case we can approximate the matrix element in the 



lowest order of {2Ec4:/E, 



J4J 



^1/4 



by 



{9\Ti\e) 



gmi 



(55) 



where we have labeled the two lowest eigenstates of the 
four-junction flux qubit as \g) and |e). In the limit of 
weak coupling £^j4 3> we find lj ^ gmi. Then the 
matrix element simplifies to {g\T^\e) — ~iaEj'mi/8Ej4^. 
In general we have £'c4/£'j4 ^ 1 and therefore transi- 
tions due to tunneling across the fourth junction are well 
suppressed compared to tunneling across the other junc- 
tions. Similarly we can calculate the tunneling through 
the third junction with the tunneling operator TJj ~ 
gi{<l>+-'t'-)/2 g^jjj^^ using the same approximation as for 
Eq. dini), we get 



{g\T,\e) « a IT3I t) 



(56) 



In conclusion we can say that the fourth junction should 
have little effect on the overall decoherence rate of the 
p-flux qubit. 



D. Tunable gap flux qubit 

In the tunable gap flux qubit the small junction is re- 
placed by a SQUID, see Fig. [Tfc). In this case the pa- 
rameter a can be tuned by changing the field ipx- In 
Hamiltonian pop this changes 



a = ao cos 



(f) 



The four tunneling operators are given by 

T, — p»('P++¥'-)/2 rp _ i(ip<,xt+Vi/2-2<p+)/2 



rp _ gi(ipcxt-¥'x/2-2ip+)/2 J. _ gi(ip+-</3_)/2 



(58) 



where Ti corresponds to the junction with phase ipi. 
Since the eigenstates are the same as for the p-flux qubit, 



it is straightforward to calculate the tunneling matrix 
elements. The decay rates for quasiparticle tunneling 
through the large junctions are the same as for the stan- 
dard persistent-current flux qubit and we therefore get 
the similar overall rates. The only difference is that 
the tunneling across the small junctions, that form the 
SQUID, have now different type of interference effect if 

iPx 7^ 0. 



E. Capacitively shunted Flux qubit 

Another version of the same qubit, the capacitively- 
shunted flux qubit— has the same Hamiltonian as the 
persistent-current flux qubit pop . The major difference 
is that it is operated in the regime a < 1/2 and the 
charging energy along the dimension of the double well 
potential is significantly reduced by introducing a shunt 
capacitance, as 3> a. This qubit has excellent coherence 
times and still preserves relatively large anharmonicity. 
We will discuss the system at the symmetry point tpcxt — 
TT and for the case where it deviates most strongly from 
a harmonic oscillator, a = 1/2. After this we are able to 
generalize this result to other parameter regimes. 

Under the conditions we specified above we can write 
the Hamiltonian as 



-ffqubit cxt 



1/2) ^E'cNl + ^E.^l, (59) 



with E^j = Ec/{1 



2a + Cs/Cj). We note that the 
Hamiltonian is similar to the one of the strongly anhar- 
monic phase qubit^. We find the eigenstates to be 



(^+U> = 



(^+lt) = 



1/4 



1/4 



(60) 



(61) 



with u' = [iEj/l&E'f^f'^. It is now simple to calculate 
the matrix elements for the tunneling operators (j32p . 



aifiit) = aiT3it) 



1 C^Jk. 

2 \iE,j 



1/6 



(62) 



and for T2 one obtains again the additional phase fac- 
tor 6*'^'='='/^ = i. The energy splitting between the 
eigenstates is given in a good approximation by 5E — 
(57) 2{12E'(j'^ Ej)^/'^ . Using this result we can rewrite the 
matrix element as (| iTij f) — \/E'fj/5E. This is a well 
known result for a harmonic oscillator and is therefore 
also valid for the complete relevant parameter regime of 
the capacitively shunted fiux qubit. The interference an- 
gles (pQj of the three junctions are the same as for the 
p-fiux qubit at the symmetry point: In the case of non- 
equilibrium quasiparticles and qubit decay the junctions 
1 and 3 have constructive {^pm = lySos = 0) and the junc- 
tion 2 has destructive interference ((^02 = tt). Vice versa 
for the dephasing. 
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As mentioned, the c-flux qubit is very similar to the 
strongly anharmonic phase qubit, with the difference that 
the latter qubit has only a single junction, which corre- 
sponds to the second junction of the c-flux qubit. As this 
junction shows destructive interference against qubit de- 
cay due to non-equilibrium quasiparticles, the strongly 
anharmonic phase qubit is better protected against this 
decoherence source. However, the matrix elements of the 
qubits are generally similar to that of harmonic oscilla- 
tors, and this means that the c-flux qubit is relatively 
well protected against quasiparticle tunneling, too. In 
table U we compare the qubit decay times between the 
capacitively shunted flux and the other discussed qubits. 



V. CONCLUSION 

In this work we analyzed decoherence in supercon- 
ducting qubits due to quasiparticle tunneling. We con- 
sidered two types of sources of quasiparticle tunneling 
relevant for low temperatures and low qubit energies: 
non-equilibrium quasiparticles and sub-gap density of 
states. Using typically observed values for their densi- 
ties we estimated the resulting qubit decay times and 
compared them with the experimentally measured ones. 
In the best case scenario, i.e., with the lowest reported 
quasiparticle densities, we showed that only the decay 
times of persistent-current flux qubits were similar to re- 
cent experiments. The multi-junction flux qubits have 
achieved recently excellent coherence times and have 
the remarkable feature of extremely large anharmonicity 
which makes it the best approximation to a real two level 
system of all superconducting qubits. However, the pre- 
sented analysis shows that such qubits are also very frag- 
ile against quasiparticle tunneling induced decoherence. 
To protect the qubits against creation of non-equilibrium 
quasiparticles, for example, a careful isolation from the 
nearby environments should be realize d^*^'^"^ . 
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APPENDIX 

A. Derivation of the decoherence rates 

Here we formulate the qubit decay rate as a function of 
electron correlation (Green's) functions of the leads. Our 
aim is to consider the response of the superconductor as 
a sum over all energy states, rather than consider the 
response of single states. By this we can model situations 



where the superconductor quasiparticle density of states 
at a given energy cannot be mapped back into certain 
quasiparticle energy in the normal state, for example, 
due to energy-level broadening effects. 

We consider a single Josephson junction qubit and 
write the time evolution of the reduced density matrix, 
when interacting with quasiparticle environment, gener- 
ally as 



Amn(i) = i{En - Em) 



f dt'a^an{t-t')Pab{t'). 

ab 



(63) 



Here the (generalized transition rate) tensor cr^J^™ in- 
cludes the effect of qubit-quasiparticle interaction. Its 
calculation is similar to the case of metallic reservoirs^ 
with a difference that when tracing out the environment 
also two annihilation (creation) operators c^^"* can also 
contract to pairs. Their contribution leads to the pair 
density p, whereas diagonal contributions, such as c^Ck, 
lead to the density of states n. 

In the leading order we obtain for the transition from 
the state |i) to |/) 



(64) 



ait-t') = {mH^{t)\f){f\H^{m) 



The trace over the initial distribution of quasiparticles 
leads to four nonvanishing contributions 



a{t) = Le»'5-*(a + & + c + rf), 
where du — SE/h and 

kla 

& = I |.) p 5] G<Jk, i' - i)G>JM - 

kla 

c = -mf){f\n) J2 ^.>/5(k, t - t')F<^ii, t' - 1), 



(65) 



klaP 

d^-{^\f^\f){f\f^\^)J2F<^{k,t'-t)F>^il,t-t'). 

klap 

(66) 

The correlation functions G and F are defined in Eq. ([3|) . 
We define their Fourier transforms as 



1 1"^ 
27r J„oo 



(67) 



The Fourier tranforms are related to the spectral densi- 
ties as in Eq. In our analysis we have omitted the 
overall phase appearing usually in _F-functions, as it is 
treated in the qubit part. In compared to Refs. [13 and 
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[3^ our definition of corresponds to (and F< to 
F<). 

I 



Note the symmetry of the equation around uj ~ 0: contri- 
bution of negative energies {lj,lu' < 0) gives similar con- 
tribution as positive energies (w, w' > 0). As all processes 
described by a correspond to given electron tunneling di- 
rection, negative energies correspond to opposite tunnel- 
ing direction of the quasiparticle. The contribution from 
b follows the same calculation and is similar to a, with 
a difference that the qubit matrix element is changed to 
K/lTtl^^p corresponds to tunneling of an electron to 

I 



where we have used the fact that for the considered qubits 
|(/|f = |(/|ft|z)| and then defined {i\f\f){f\f\i) = 



J 



We 



leading to equation (fT51) . This expression is similar to 
equation (9) in Ref. |34| . giving also insight to the inter- 
pretation of the different terms. The main term that 
is missing, when compared to Ref. is the coherent 
Josephson term (sint^-term), as it does not contribute 
in the process considered here. (This term is treated 
exactly in the qubit Hamiltonian) . Such coherent terms 
would contribute through the principal value integration. 



Provided by the assumption of constant tunneling am- 
plitude t we can now write the contributions to the tran- 
sition rate coming from the different terms as 



(68) 
TrS{uj) — 



(69) 

I 

the opposite direction. 

The contribution from c and d exists only if a factor 
{i\f\f){f\f\i) = (z|ft|/)((/|ft|i))* is finite. It means 
that it exists if two electron tunneling directions can lead 
to the same final state, in a single quasiparticle tunneling 
process. Then the processes interfere. The contribution 
from c is similar as before, but replaces n functions by p 
fimctions, 




I 

-e''^°\{f\f\i)\'^. Noticing, that the contribution from d 
is equal to contribution from c, except with an opposite 
phase factor e"*'''", one obtains that the principal value 
part gives no contribution. Therefore one gets 



I 

The dissipative terms are similar. Especially, the famous 
cos (y9-term corresponds to interference between electron 
and hole-like tunneling in our calculation, described by 
terms c and d. A new type of measurement of this term 
with a superconducting charge qubit has been considered 
in Ref. [11 



:^|(/|f limRe<| e^''"*' / dt'e~^'~*'>l duj I doj'n(ij)n(uj')f+(u})f-(uj')e''^'^'-'^'^'' 



We have defined the normalized density of states n as in Eq. ([S|). Using the relation limg^o jl^ g^i^t+st 
iP{l/{uj)), one obtains 

/"OO />oo 

a-^ D^TTt^\{f\f\i)\^ / duj duj'f+{uj)r{uj')n{uj)n{uj')5{uj-uj' -5uj). 



OO J — OO 



ft2 



\{f\m\'Re\e^^o 



duj 



du'J+{u:)r{^')p{^)p{^') 



5{lj ~ lj' ~ lj fi) — 



c + d^ \{f\T\i)\^ COS / du duj'f+{u;),riuj')p{Lj)piu;')Siu;-uj' -Uf,). (71) 

can write now 

a + b + c + d^ —:^\{f\f\i)f duj dLj'f+{io)f~{Lu')[n{uj)n{uj)+cosipop{io)p{uj')]S{io^uj' -ujf,), (72) 

" J — OO J — OO 
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B. Coherent terms: Parameter renormalization 

Taking into account the tunnehng Haniihonian by 
perturbation theory includes not only incoherent pro- 
cesses but also coherent terms leading to renormaliza- 
tion of the qubit parameters. Their contribution is usu- 
ally small but can depend, for example, on the non- 
equilibrium quasiparticle density. This has been studied 
in detail in Ref. [l2|- 

The Josephson coupling i?j in Eq. (2) is calculated as 
the expectation value over the quasiparticle distribution 
of the operator 



El 

2 



1 



{N + 2\Ht Ht\N) 



qp 



(73) 



Here the quasiparticle (BCS) Hamiltonian has no depen- 
dence on the charge number N. Using fermi distribu- 
tions for the excitation occupation probabilities one ob- 
tains the famous Ambegaokar-Baratoff relation for Ej. 
This corresponds to the usual choice of the Josephson 
coupling. In addition, there exists a similar second-order 
contribution describing reactive behaviour of quasiparti- 
cles (usually referred to as the capacitance renomaliza- 
tion) 



He = \N){N\Ht^Ht\N){N\ 



(74) 



To calculate this one can introduce a cut-off function 
D{E) = 1/[1 -I- (E/Eco)"^], where Eco is larger than any 
of the relevant energy scales in the system introduced to 
avoid divergence of the energy corrections. At this point 
this term produces a constant energy shift for all N and 
therefore does not bring anything new on the system. 

Insted of including these terms into the time- 
dependent problem, it is easier to estimate their effect by 
considering energy-level changes using time-independent 
perturbation theory. Using the second-order theory 
again, but now for the modified tunneling Hamiltonian 
i?T , one obtains a correction to the energy of the state i 



SH-\ = Ej{i\ cos.^|z) 



+ ( mT 



1 



Ei -^qp -^qubit 



(75) 



where Ei is the energy of the qubit states At this 
point we see that, for example, the Ambegaokar-Baratoff 
result is a good approximation as long as the approxima- 
tive Josephson term _Ej cos ip and the corresponding term 
coming from the second term on the right-hand side of 
Eq. ([75]) almost cancel each other. It occurs if the tran- 
sition probabilities to virtual states \v) {v ^ i) are small, 
or if the energies of the virtual states stay small i?„ ^ A. 
In the opposite case extra anharmonicity (non-constant 
energy-level spacing) could emerge due to this correction, 
as the effective tunneling coupling can become different 
for different states. 

C. Decoherence due to higher-order effects: 
Andreev tunnehng 

The decoherence processes / and // correspond to 
leading order (incoherent) tunneling effects in the system. 
The most important process in higher orders is usually 
the Andreev tunneling. This involves, for example, tun- 
neling of two electrons from one side, with leaving two 
excitation behind, to form a Cooper pair on the other 
side. Such processes dominate the sub-gap conductance 
of normal metal-insulator-superconductor junctions for 
ideal BCS state on the superconductor side. However, 
for the case of perfect SIS junctions no such process exist 
until the energy 2A is somehow provided for creation of 
two excitations. The process appears as a step-like be- 
haviour in the current- voltage characteristics of single J Js 
nearby eV = A. In experiments the height of the step 
is usually considerably higher than obtained by theory. 
The reasons for this are still unclear—. 

For the case of superconducting qubits no energy is 
available to create two excitations, unless sub-gap den- 
sity exists at one side of the junction. Through creation 
of excitations into subgap region with creating a Cooper- 
pair on the other side of the junction one obtains a pos- 
sible contribution to the decoherence. However, a simple 
analysis indicates that the rates for such processes are 
proportional to (r/A)^/A^, where ^ 1 is the effective 
number of parallel tunneling channels. It results that 
contribution from this process should be much smaller 
than of the process // considered in this paper. We con- 
clude that the contribution from higher-order tunneling 
effects to decoherence in superconducting qubits should 
stay small. 
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